16s分析之进化树+差异分析（一）
之前我便说过，16s出图主要是用R，自这两年来Y叔的ggtree出来，好评不断，引用不断攀升
这里当然借用Y树的包来做我们的进化树，当然Y叔公众号讲的更加深入和详细，有需要请关注：biobabble
这里Y树曾经详细介绍过ggtree，里面还顺带介绍了ggplot，有兴趣的亲查看（反正我收获颇丰）：
这里我参考了宏基因组公众号的代码：（参考别人的嘛，也不难为情，工具这种事情，总是不多看别人怎么做，自己学着搞）
这是刘永鑫的博客：http://blog.sciencenet.cn/blog-3334560-1079159.html
现在我们开始：
这里我的R包很老了，Y叔的包还是很新的，不多说，只能更新R包了
##################首先我更新了一下R，已经不能安装ggtree###################
#由于R太老了，所以我跟新了R
#安装更新工具
install.packages("installr")
library(installr)
updateR()
#安装ggtree
source("https://bioconductor.org/biocLite.R")
biocLite("ggtree")
#更新的R包有好多依赖包没有，手动安装，还有好多包需要升级，慢慢弄吧
source("https://bioconductor.org/biocLite.R")
biocLite(c("treeio"))
install.packages("tidyr")
##################首先我更新了一下R，已经不能安装ggtree###################
下面开始我们的分析，在画进化树之前首先我们的前处理在qiime中：
# 选择OTU表中丰度大于0.1%的OTU，足够我们分析的了
filter_otus_from_otu_table.py--min_count_fraction 0.001 -i otu_table.biom -o hg_tree/otu_table_k1.biom
# 获得对应的fasta序列rep_seqs.fa
filter_fasta.py-f rep_seqs.fa -o hg_tree/rep_seqs_k1.fa -b hg_tree/otu_table_k1.biom
# 统计序列数量，正则表达式，如果看不懂，百度搜索，一大堆，这里我简单说一下，单引号中即为所搜索的字符，-c就是计算找到 '搜寻字符串' 的次数，输出含有>的次数
grep -c '>'hg_tree/rep_seqs_k1.fa
这里我们安装clustalo
#qiime已经集成了安装多序列比对的依赖包了，我下载的版本：clustalo-1.2.4-Ubuntu-x86_64，使用qiime1.9.2亲测可直接安装此包，无需依赖
sudo apt-get installclustalo
# clusto多序列比对
clustalo -ihg_tree/rep_seqs_k1.fa -o hg_tree/rep_seqs_k1_clus.fa --full --force--threads=5
# 使用qiime的脚本调用fastree建树
make_phylogeny.py-i hg_tree/rep_seqs_k1_clus.fa -o hg_tree/rep_seqs_k1.tree
# 格式转换，这里我们要提一下，为什么要做格式转换，这就是我们使用qiime的make_phylogeny.py命令做的树文件，首先它是一个无根树，其次，最重要的是，每个OTU名称都被单引号夹起来了，当然我们不需要这种形式，所以下面开始去除
到这里，可能就有很多人都有满肚子疑问了，tree的文件到底是一种什么格式，我们这种格式得给的通俗的解释啊！别着急，下面我来做这个工作：
Newick格式，就是我们使用qiime做进化树的格式，也是最为古老的格式，是一种使用圆括号和逗号代表具有节点长度理论图形树的一种方法。
library("ggtree")
library("ggplot2")
#只命名叶节点
tree1 <-"((A, B), C,D);"
x1 =read.tree(text = tree1)
ggtree(x1)+geom_tiplab()
#没有任何一个节点被命名
tree2 <-"((,),,);"
x2 =read.tree(text = tree2)
ggtree(x2)+geom_tiplab()
#所有节点均命名
tree3 <-"((A, B)E, C,D)F;"
x3 =read.tree(text = tree3)
#这里我显示叶节点并调整了叶节点标签,但是如何显示节点标签，我就不会
ggtree(x3)+geom_tiplab(offset=1,size=2,align=1)
#
#除了根节点，所有节点均命名
tree3 <-"((A, B)E, C,D);"
x3 =read.tree(text = tree3)
#这里我显示叶节点并调整了叶节点标
ggtree(x3)+geom_tiplab(offset=1,size=2,align=1)
#这里给加上了节点之间的距离
tree3<-"(A:0.1,B:0.2,(C:0.3,D:0.4)E:0.5)F;"
x3 =read.tree(text = tree3)
#这里我显示叶节点并调整了叶节点标签,但是如何显示节点标签，我就不会
ggtree(x3)+geom_tiplab()+xlim(NA, 1) + ylim(NA, 5)#offset=1,size=2,
ggtree(x3)+geom_tiplab()+xlim(0, 1) + ylim(0, 5)
我们在上面简单了解了一下树的知识，下面我们来做一颗树
# 这次我使用番茄otu
setwd("E:/Shared_Folder/HG_usearch_HG/ggtree_hg")
# 读入分析相关文件
# 读取树文件
tree3 <-read.tree("rep_seqs_k1.tree")
# 读取树物种注释信息
tax <-read.table("rep_seqs_k1g.txt",row.names=1)
head(tax)
tax$id=rownames(tax)
head(tax)
# 物种注释等级标签，六级
colnames(tax) =c("kingdom","phylum","class","order","family","genus","id")
## 给每个OTU按门分类分组
groupInfo <-split(row.names(tax), tax$phylum)
tree3 <-groupOTU(tree3, groupInfo)
我们来出一张默认树图看看，可以看到树建立在活动窗口，当前设备，这是一颗无根树
ggtree(tree3)
#使用branch.length= "none"，叶对齐，默认不是对齐的
ggtree(tree3,branch.length = "none")
#调整树的结构ladderize= F
ggtree(tree3,branch.length = "none",ladderize = F)
#改变树的形式，按照跟组上色
ggtree(tree3,layout="fan", branch.length = "none",ladderize = FALSE,aes(color=group))
#调整树的开合角度和按照一定角度旋转
p<-ggtree(tree3,layout="fan", branch.length = "none",ladderize = FALSE,aes(color=group))
open_tree(p,angle=100)%>% rotate_tree(-45)
#我们是在没有必要写这么多注释，这里我们意思一下
p<-ggtree(tree3,layout="fan", branch.length = "none",ladderize = FALSE,aes(color=group)) %<+% tax1+
geom_tiplab2(size=1.5)+ theme(legend.position= "right") +
geom_tiplab2(aes(label=phylum,color=group),size=2, offset=6)
#查看节点，为什么查看这个，因为我们要按照节点给他高亮
p+geom_label2(aes(subset=!isTip,label=node),size=2)
##加上阴影
p1<-p+
geom_hilight(node=323, fill="gray",alpha=0.5) +
geom_hilight(node=320, fill="pink",alpha=0.5) +
geom_hilight(node=313,fill="beige", alpha=0.5) +
geom_hilight(node=298,fill="yellow", alpha=0.5)+
geom_hilight(node=255, fill="blue",alpha=0.5) +
geom_hilight(node=237, fill="red",alpha=0.5) +
geom_hilight(node=194,fill="green", alpha=0.5)+
geom_hilight(node=269,fill="lightblue", alpha=0.5)
p1
今天的推送就到这这里，但是没有完，树的这些东西还有好多，我正在准备，下次发
这两天的确是发文不如之前迅猛了，不是疲软，还没有到那个年纪，只是在做一些其他的分析而已
最后还是那句话：
学习永无止境，分享永不停歇！